Cite this Article
Xie Qing, Wang Lei, Li Jiangxu, Li Ronghan, Chen Xing-Qiu. General principles to high-throughput constructing two-dimensional carbon allotropes. Chinese Physics B, 2020, 29(3): 037306  
Permissions
General principles to high-throughput constructing two-dimensional carbon allotropes
Xie Qing1, 2, Wang Lei1, 3, Li Jiangxu1, 3, Li Ronghan1, 3, Chen Xing-Qiu1, †
Shenyang National Laboratory for Materials Science, Institute of Metal Research, Chinese Academy of Sciences, Shenyang 110016, China
University of Chinese Academy of Sciences, Beijing 100049, China
School of Materials Science and Engineering, University of Science and Technology of China, Shenyang 110016, China

 

† Corresponding author. E-mail: xingqiu.chen@imr.ac.cn

Project supported by the National Science Fund for Distinguished Young Scholars, China (Grant No. 51725103) and the National Natural Science Foundation of China (Grant No. 51671193). All calculations have been performed on the high-performance computational cluster in the Shenyang National University Science and Technology Park.

Abstract

We propose general principles to construct two-dimensional (2D) single-atom-thick carbon allotropes. They can be viewed as the generalization of patterning Stone–Walse defects (SWDs) by manipulating bond rotation and of patterning inverse SWDs by adding (or removing) carbon pairs on the pristine graphene, respectively. With these principles, numerous 2D allotropes of carbon can be systematically constructed. Using 20 constructed 2D allotropes as prototypical and benchmark examples, besides nicely reproducing all well-known ones, such as pentaheptites, T-graphene, OPGs, etc, we still discover 13 new allotropes. Their structural, thermodynamic, dynamical, and electronic properties are calculated by means of first-principles calculations. All these allotropes are metastable in energy compared with that of graphene and, except for OPG-A and C3-10-H allotropes, the other phonon spectra of 18 selected allotropes are dynamically stable. In particular, the proposed C3-11 allotrope is energetically favorable than graphene when the temperature is increased up to 1043 K according to the derived free energies. The electronic band structures demonstrate that (i) the C3-8 allotrope is a semiconductor with an indirect DFT band gap of 1.04 eV, (ii) another unusual allotrope is C3-12 which exhibits a highly flat band just crossing the Fermi level, (iii) four allotropes are Dirac semimetals with the appearance of Dirac cones at the Fermi level in the lattices without hexagonal symmetry, and (vi) without the spin–orbit coupling (SOC) effect, the hexagonal C3-11 allotrope exhibits two Dirac cones at K and K′ points in its Brillouin zone in similarity with graphene.

PACS: ;73.20.At;;63.22.Np;;73.22.Pr;
Keyword:;graphene allotropes;;phonon spectra;;semimetal;;dirac Fermions;;modeling;
1. Introduction

The experimental realization of graphene has stimulated tremendous research activities on this fascinating two-dimensional (2D) one-atom-thick material[1–13] because of its revolutionary impacts to future industry in many fields, such as energy storages,[15–18] solar cells[19,20] nanoelectronics,[21,22] as well as superconductor.[23–25] Nevertheless, these various applications require versatile and controllable properties upon modifications of the graphene.

Within this context, many two-dimensional single-atom-thick carbon allotropes, called graphene allotropes (GAs), have been proposed. Even before the graphene was confirmed experimentally, one already paid much attention to 2D carbon allotropes. As early as in 1987, some graphynes were theoretically suggested[26] and they were recently demonstrated in details to exhibit the direction-dependent Dirac cones through first-principles calculations.[27,28] In 1994, some other 2D carbon networks were designed for the conjugated-circuit computations.[31] In 1996, one of those proposed 2D networks,[31] later called pentaheptite,[32] composed of pentagons and heptagons, was theoretically claimed to be a metal but with covalent carbon bonds. Just one year later, a 2D net-C structure, composed of the 4y+6+8 membered rings of carbon, was proposed,[33] while in 2000 the Haeckelite structures containing 5+6+7 membered rings of carbons were further presented.[34] Most recently, it was shown that the family of the 5+6+7 membered rings of carbon harbors distorted Dirac cones.[35]

Since graphene was confirmed experimentally,[1,2] the search and study of new GAs was spurred again. Through various nano-engineering defects on the pristine graphene, a class of memberanic carbon allotropes, named dimerite, were conceived.[36,38] Subsequently, based on the defect patterns some 2D carbon allotropes have been proposed,[39] such as octite-SC (5+6+8 membered rings),[40] T-graphene (4 + 8 membered rings),[41] two OPGs (5+8 membered rings),[44] tetrahexcarbon,[92] C10,[93] and net-W (4+6+8 membered rings).[45] In addition, using high-throughput structural searching codes (e.g., CALYPSO[46] and USPEX[47]), numerous GAs (e.g., S- and T-graphene) have also been constructed.[35,48] Recently, a new allotrope with only pentagons, called penta-graphene,[49,50] was proposed. Although many different GAs have been reported, there is still lack of a simple and systematic route to effectively seek structural candidates.

Physically, a fascinating property of graphene is the presence of Dirac cones at the Fermi level. Therefore, the low-energy excitations in graphene are massless Dirac fermions. This peculiar quasiparticle is responsible for many exotic phenomena observed in graphene, for examples, fractional quantum Hall effect[3,14] and ultrahigh carrier mobility.[51] In connection with similar massless Dirac fermions in graphene, it is another frontier of topological Dirac semimetals most in 3D bulk materials, such as Na3Bi and Cd3As2 families,[52–54] and of topological Dirac nodal line semimetals featured by a closed curve by band crossing in analogy with a closed ring formed by Dirac cones, found both in pure metals[55,58,74] and a series of interesting compounds.[56–70] They also exhibit 3D Dirac cones around the Fermi level, inducing exotic properties.[71–78] However, very few 2D Dirac materials are known (see Ref. [79] for a review). The reason is that Dirac cones are topologically protected and robust against perturbations, thus, cannot be obtained by fine tuning the Hamiltonian parameters. This makes discovery of 2D Dirac materials a nontrivial task.

In this paper, we propose general schemes of altering the carbon structure on the pristine graphene. The first one is to implement bond rotation operation, which is a generalization of patterning Stone–Walse defects (SWDs) on the pristine graphene. The second scheme is adding (removing) carbon pairs. This can be viewed as a generalization of patterning inverse SWDs (iSWDs) on the pristine graphene. Based on these two general schemes, we construct numerous 2D carbon networks, including many well-known ones such as pentaheptites,[32] pentahexoctite,[42] T-graphene,[41] OPGs,[44] and so on.[29,43] Moreover, we report many new structures with combinations of other types of membered rings. For these new structures, the formation energies are found to be higher than that of graphene within 1.2 eV per carbon atom. Therefore, they are thermodynamically metastable. However, most of these structures are dynamically stable, due to the absence of imaginary modes in their phonon spectra. Two structures are dynamically metastable. The electronic band structures show that their exhibit various conducting behaviors ranging from metallic, semimetallic to semiconducting. Interestingly, there are four new GAs with massless Fermi fermions due to the presence of Dirac cones at the Fermi level. The C3-11 structure with hexagonal symmetry exhibits an isotropic linear dispersion near K(K′) point, which mimics the case in graphene. The other three new GAs without hexagonal symmetry have anisotropic linear low-energy excitations. These results dispute the long-standing belief that the presence of Dirac cones in graphene is related to its hexagonal symmetry. Our results show that the 2D carbon system is a fertile ground for searching for Dirac materials.

2. Two constructing schemes

It is well-known that there are two fundamental defects on graphene: the SWD[37] and the iSWD.[38] The formation of a SWD can be understood by rotating the interior carbon–carbon bond in a diamond unit grouped by four adjacent hexagons (Fig. 1(a)). An iSWD can be formed by adding two extra carbon atoms onto the centering hexagon of the basic unit composed by three line-up hexagons (Fig. 1(b)). As early as 1990s, patterning SWD on pristine graphene was demonstrated to be an efficient method to construct 2D planar networks.[36,38] These two defects correspond to two basic operations: the bond rotation and the addition of extra atoms. Note that in a SWD each bond rotation operation transforms four hexagons into two pentagons and two heptagons. But, in an iSWD the addition of the extra two atoms transforms three line-up hexagons into two pentagons and two heptagons.

Fig. 1. (a) A SWD is formed by rotating the interior bond of a diamond unit. (b) An iSWD is formed by adding two extra atoms onto the parallel bond pair. By adding two extra atoms to the 2nd and 1st bond pairs, a 4-atom ring in (c) and a 3-atom ring in (d) are formed, respectively.

Here, we propose the first scheme by allowing performing bond rotation on four adjacent any membered carbon rings. The effect of this operation is to extend the two adjacent rings by one atom (i.e., from N1− and N2−membered to N1 + 1 and N2 + 1 -membered rings) while shrink the non-adjacent two rings by one (i.e., from M1 - and M2 -membered to M1 — 1 - and M2 — 1-membered rings). A SWD is a special case where the bond rotation is performed on four adjacent hexagons (i.e., N1 = N2 = M1 = M2 = 6).

There are six carbon–carbon bonds in a hexagon on the pristine graphene. These six bonds form three types of bond pairs: 1) two bonds (of 120°) that share a carbon atom, 2) two bonds (of 60°) that are connected by another carbon bond, 3) two bonds (of 0°) that are parallel to each other. An iSWD is formed when two extra carbon atoms are added to the third kind of bond pair (see Fig. 1(b)). We propose a second scheme by allowing adding/removing extra carbon atoms to/from different locations. Two additional examples are given in Figs. 1(c) and 1(d), where two extra carbon atoms are added to the 2nd and 1st kind bond pairs, respectively. By exploiting these two general schemes, numerous 2D carbon networks can be systematically constructed. In brief, these two schemes can be described as follows:

The first scheme: A SWD is formed by performing bond rotation on four six-membered rings. The first scheme generalizes this to allow performing bond rotation on four any-membered rings.

The second scheme: An iSWD is formed by adding two extra atoms to a specific location on the hexagon. The second scheme is to allow adding extra atoms to many different locations.

3. Methods

We perform first-principles calculations based on density functional theory with generalized gradient approximation (GGA) in the form of Perdew–Burke–Ernzerhof[81] exchange–correlation functional. All the calculations are performed using Vienna ab initio simulation package (VASP)[82] with projector augmented wave method (PAW).[83] The vacuum slabs of 15Å are chosen to avoid interactions between adjacent carbon layers. A self-consistent field method (tolerance eV/atom is employed in conjunction with the plane-wave basis sets of cutoff energy of 600 eV to obtain very accurate results. Geometry optimization is implemented until the remanent Hellmann–Feynman forces on each ion are less than 0.001eV/. We use a Γ-centered k point mesh to sample the Brillouin zone (BZ) for the structures with hexagonal unit cells. For other unit cells, we use the Γ-centered Monkhorst–Pack k point mesh.[84] Phonon spectra are obtained using the finite displacement method. A 4 × 4 × 1 supercell is used for the calculation of the phonon spectra to make sure that the force constants are sufficiently collected. The free energies as a function of temperature are calculated by using phonon dispersion within the quasiharmonic approximation. All the phonon spectra, vibrational density of states (DOS), and free energy are calculated using Phonopy package.[85]

4. Results and discussion
4.1. Construction of GAs

The 4+8 membered rings. On pentaheptites, we apply bond rotation upon the interior bond of an iSWD unit whose boundaries are defined by the yellow bonds in Figs. 2(b) and 2(c). Then the two pentagons are both declined to four membered rings (C4), whereas the two remaining heptagons are extended to octagons (C8). Assuming that bond rotations are simultaneously performed on all iSWDs patterns on pentaheptite-I or pentaheptite-II, the so-called T-graphene[41] (Fig. 3(d)) can be constructed, which features combinations of C4 and C8.

Fig. 2. The optimized structures of graphene allotropes. The yellow bonds in (a)–(c), (f), (g), (i), (j) define the boundary of building blocks. The red bonds are to be rotated in order to construct new structures (see the context). The red carbon atoms in (l)–(n) denote the extra ones. The red carbon atoms in (e), (k), (o)–(t) mark the differences with triangle units. The dashed blue lines mark the primitive cells.
Fig. 3. (a) An armchair ribbon with line-up diamonds. (b) A zigzag ribbon with zigzag arranged diamonds. The dashed cycles denote diamonds. Two neighbor diamonds share one hexagon. By rotating the interior red bonds of all diamonds, the sharing hexagons become octagons (8) and the non-sharing hexagons become pentagons (5). Finally, two strips of pentagons and octagons are constructed.

The 3+9 membered rings. There are no common bonds between any two C4 rings on T-graphene, one can not obtain a 3+9 membered ring straightforwardly. However, through clever design, it is still possible. One realization is presented in Fig. 2(d). By rotating the red bonds on T-graphene, a structure with 3+9 membered rings can be constructed (see Fig. 2(e)).[31,43] Throughout this transformation, the final result is that each C8 is transferred to a C9 by increasing one atom, while all C4 are reduced to C3 by subtracting one atom. Note that a C8 ring is effected by nearby three bond rotations, two of which each increases the C8 ring by one atom. The other one decreases the C8 ring by one atom. Therefore, the final result is to extend C8 to a C9. Certainly, one may proceed along this route to continue creating other 2D allotropes by reducing C3 to acetylenic C2 linkages (–C=C–). Then numerous graphynes also can be constructed.[27–30]

The 5+8 membered rings. Geometrically, graphene can be divided by other types of diamonds combinations. For instance, we consider the cases where two neighboring diamonds commonly share a hexagon as illustrated in Figs. 3(a) (an armchair ribbon with line-up arrangement of diamonds) and 3(b) (a zigzag ribbon with zigzag arrangement of diamonds), respectively. Using these arrangements, after rotating the interior red bonds of every diamond, every commonly shared hexagon is extended to an octagon and non-sharing hexagons are shrink to pentagons. We obtain two different ribbons (one with armchair edge and one with zigzag edge) composed of pentagons and octagons. By repeating these ribbons along vertical directions, a class of allotropes, called OPGs, can be constructed. From Fig. 2(a), one unique OPG (OPG-L[44] in Fig. 2(f)) can be created. By translating the basic ribbon in Fig. 2(b), two geometrically inequivalent OPGs (OPG-Z[44] in Fig. 2(g) and OPG-A in Fig. 2(h)) can be constructed. One can create other OPGs by simply employing more ribbons as the repeating units.

The 4+10 and 3+12 membered rings. Based on OPG-L and OPG-Z, 4+10 membered rings can be constructed by rotating the common bonds of two adjacent pentagons. Figures 2(i) and 2(j) present two simplest C4-10-I and C4-10-II structures obtained by rotating the red bonds in Figs. 2(f) and 2(g). We note that from OPG-L, many geometrical inequivalent C4-10 allotropes can be constructed depending on which common bonds of the neighbor pentagons are to be rotated. This is because in OPG-L, the pentagons are lined up forming a pentagon ribbon. There are many ways to choose the rotated common bonds to extend the neighbor C8 to C10. While in OPG-Z, the pentagons are zigzag arranged to form a pentagon ribbon, and in OPG-A, four pentagons are surrounded by six octagons. Furthermore, rotating the common bond of two neighboring C4 rings in C4-10-I (the red bonds in Fig. 2(i)) and C4-10-II (the red bonds in Fig. 2(j)) yields a unique hexagonal C3-12 in Fig. 2(k), which is composed of 3+12 membered rings. Up to now, we find that the SWD can be used to change the structure from orthorhombic lattice to hexagonal lattice, such as from T-graphene to C3-9-H and from C4-10-II to C3-12, which is closely associated with the hexagonal graphene and four bonds used in the SWD method.

The 4+7, 4+8, and 4+12 membered rings. It is heuristic to compare the geometries among graphene, T-graphene, C4-10-I, and C4-10-II. By removing two common carbon atoms of adjacent C4 rings in the C4-10-I or C4-10-II structures, two adjacent C4 rings decrease to one. This eventually yields T-graphene. On the other hand, the pristine graphene can be recovered by removing two carbon atoms of each C4 ring on T-graphene. The above observations demonstrate the usefulness of our second scheme in connecting the known structures. We now show that it is also very useful in constructing new structures. For instance, by adding extra carbon pairs (the red atoms in Fig. 2(l)) on graphene as illustrated in Fig. 1(c), a 2D allotrope (called C4-7) composed of C4 and C7 rings is constructed in Fig. 2(l). By removing carbon pairs in every two C4 cycles in T-graphene, the net-W and net-C structures containing 4+6+8 membered rings can be constructed, as reported in Ref. [45]. On the net-W allotrope, adding carbon pairs (the red atoms in Fig. 2(m)) results in an orthotropic C4-8-R structure composed of C4 and C8 in Fig. 2(m). Similarly, another C4-12 structure in Fig. 2(n) containing C4 and C12 rings can easily be constructed from C4-10-I (the red atoms mark the differences).

The 3+7, 3+8, 3+9, 3+10, and 3+11 membered rings. We compare the structures of pristine graphene (Fig. 2(a)), C3-9-H (Fig. 2(e)), and C3-12 (Fig. 2(i)). They all have hexagonal symmetry. Furthermore, they exhibit close geometrical relevance. The C3-12 structure can be constructed via replacing all carbon atoms (i.e., A-sublattice and B-sublattice) by triangle units following the scheme in Fig. 1(d). The C3-9-H structure emerges by replacing half carbon atoms (i.e., A-sublattice or B-sublattice) by triangular units. Inspired by this observation, we can choose different carbon atoms to be replaced by triangle units. For instance, by replacing every common atom of three adjacent hexagons with triangular units (this transforms the three adjacent hexagons C6 to heptagons C7) on graphene, numerous C3-7 allotropes can be constructed. Different C3-7 structures correspond to different patterns of triangular diamonds (a three-adjacent-hexagon unit is called a triangular diamond ) on graphene. Figure 2(o) shows the simplest C3-7 allotrope with the smallest unit cell. Similarly, we can also construct C3-8 in Fig. 2(p), C3-9-R in Fig. 2(q), C3-10-H in Fig. 2(r), C3-10-R in Fig. 2(s), and C3-11 in Fig. 2(t). They can be obtained by replacing partial atoms on graphene by triangle units. In other words, they can be obtained by replacing triangle units in C3-12 with carbon atoms. In short, they are intermediate cases between graphene and C3-12.

Up till now, we are only looking for new GAs with two types of rings (i.e., 5+7, 4+8, 3+9, 5+8, 4+7, 4+10, 3 + 11, 3 + 12). If more types of rings (i.e., 5 + 6 + 7,4+6+8, etc.) are further taken into consideration, the family of 2D GAs will be significantly enlarged. Following our two schemes, it can be easily done.

4.2. Properties

Energetic stabilities. We perform first-principles calculations on some selected structures considered here. Our main results are summarized in Table 1. It can be seen that the optimized lattice constants and formation energies of graphene, pentaheptites (I and II), T-graphene, OPG-L, OPG-Z, C4-7, and C3-9 are well consistent with previous reported values.[29,41,43–45] There are two allotropes (T-graphene and C4-7) crystallizing in square unit cells and seven in hexagonal unit cells, as illustrated in Fig. 3. In particular, graphene has the smallest hexagonal unit cell with two carbon atoms. The pentaheptite-II and C3-11 allotropes contain 16 carbon atoms in their primitive unit cells. We notice that those structures with the same types of carbon membered rings but different patterns have relatively comparable formation energies. For instance, the energy difference at the ground state between pentaheptite-I and II phases is only 0.01 eV per carbon atom. Similar trend can been seen for OPG-L and OPG-Z (0.04 eV/carbon), C4-10-I and C4-10-II (0.03 eV/carbon), and C3-10-H and C3-10-R (0.03 eV/carbon). In particular, the formation energy of the OPG-A is much higher than that of both OPG-L and OPG-Z mainly due to that the density of the pentagons in the OPG-A is more accumulated (corresponding to a higher stain energy). The calculations also demonstrate that, as expected, graphene (ideal hexagons) exhibits the lowest formation energy. In fact, we also see a general way that with introducing the distorted 5+7, 4+8, and 3+9 membered rings in the ideal hexagons, the resultant structures become energetically less and less stable. This fact is apparently because the deformation energy costing for the introduction of the 5+7 rings is lower than that of the 4+8 rings with respect to the ideal hexagons and similarly, the energy cost for the 3+9 rings is higher than them. It is the key origin as to why the formation energy goes down in the sequence from graphene (−9.23 eV), pentaheptite-I (−8.99 eV), T-graphene (−8.71 eV) to C3-9 (−8.54 eV) allotropes, as evidenced in Table 1. Interestingly, following this sequence we also see that their densities become slighter and slighter as their enthalpies go less and less stable.

Table 1.

The 2D space group (2sg), Pearson symbol (Ps), optimized DFT lattice constants a, b, angle γ, atom density per area ρ (in atom/Å2), and the formation energy E (in eV/atom) of the selected 2D graphene allotropes. Z counts the number of atoms in the unit cells. E g is the band gap for semimetals and semiconductors. Y and N in the dynamic column mean dynamically stable and unstable, respectively. is the in-plane Debye temperature.

.

Vibrational stabilities. In order to study the dynamic stabilities of these selected graphene allotropes, we calculate the corresponding phonon dispersions and vibrational density of state at ambient pressure in Fig. 4. It can be seen that the phonon spectra of the graphene, T-graphene, OPG-L, OPG-Z, C4-7, and C3-9-H allotropes are in nice agreement with previous studies.[29,41,43–45] Except for OPG-A and C3-10-R structures, all other structures are dynamically stable because all their branches in the whole BZ region exhibit positive frequencies without any imaginary modes. These results imply that these allotropes are quenchable at ambient pressure conditions, once they are synthesized. It is emphasized that, due to their quite soft ZA mode, pentaheptite-II, C4-12, C3-8, C3-10-R, and C3-11 are probably unstable in the long wavelength limit,[87] which indicates that these five allotropes are sensitive to both tensile and compressive strains. In this sense, it is crucial to choose the right substrate to grow these five allotropes. Furthermore, given the importance of thermodynamical properties of 2D materials under finite temperatures, it is necessary to evaluate the free energy as a function of temperature. The free energy of crystal at a given volume is the sum of the static total energy and the vibrational free energy. Here, the static total energy at a given temperature can be obtained by fitting the curve to the third order Birch–Murnaghan equation of state. The vibrational free energy can be obtained from phonon spectra based on quasi harmonic approximation (QHA) as , where is the number of degrees of freedom in the unit cell and is Avogadro′s constant. Finally, we derive the free energy as a function of temperature (see Fig. 5). In the temperature range from 0 K to 1000 K, the graphene has the lowest free energies as compared with other allotropes and the mechanical unstable C3-9-R is the least stable thermodynamically, as well. Because using quasi-harmonic approximation above the Debye temperature is questionable, we calculate the in-plane Debye temperature of all structures by a modified 3D Debye temperature code (MechElastic[86]), as shown in Table 1.

Fig. 4. Phonon spectra and vibrational density of states (VDOS) of graphene allotropes. Two of them ((h) OPG-A and (r) C3-10-H) show imaginary modes.
Fig. 5. The free energy as a function of temperature for GAs. The data for (h) OPG-A and (r) C3-10-H are missing because they are not dynamically stable.

The Debye temperatures of graphene (2149 K) and Mos2 (610 K) are in good agreement with the values reported in Refs. [95] (2100 K) and [94] (600 K). The effective thickness of all the allotropes is 3.4 Å, which is the distance between two graphite layers. And the bulk and shear moduli of these allotropes are calculated by the isotropic modeling reported in Ref. [80]. Our calculations reveal that the C-11, C3-10-R, and pentaheptite-II allotropes become energetically more favorable than graphene at high temperatures above 1060 K, 1100 K, and 1130 K, respectively. It may suggest the high-temperature phase transformation of the graphene.

Electronic properties. The calculated electronic band structures and the densities of states for all structures are presented in Fig. 6. These allotropes exhibit various electronic behaviors ranging from metallic, semi-metallic to semiconducting. The band gaps are also summarized in Table 1. The obtained electronic band structures for graphene, pentaheptite-I, T-graphene, OPG-L, OPG-Z, C3-9, C3-12, and C4-7 allotropes are in accord with previous results.[34,41,43–45] Three notable types can be obtained as follows.

Fig. 6. Band structures and density of states of graphene allotropes. We have found (p) one semiconductor and (i), (j), (q), (t) four new semimetals. The rest structures are metallic.

1) Indirect-gap semiconductor of C3-8. Figure 6(p) compiles the DFT electronic band structure of the C3-8 structure. It exhibits an indirect band gap with the minimum band gap between high-symmetry points K and M as large as 1.04 eV. This value is an order of magnitude larger than the DFT-predicted gap of the octite SC, which can also be created by pattern defects on graphene.[40] Because the accurate band gap is important to the application of a semiconductor, we further calculate the band gap of C3-8 by employing the nonlocal HSE06 hybrid functional. As shown in Fig. 6(p), the blue bands along K−M high symmetry line are the HSE06 bands. We obtain a 1.50 eV band gap of C3-8, which is much larger (by 0.83 eV) than that of anisotropic-cyclicgraphene.[91]

2) Dirac semimetals of OPG-Z, C4-10-I, C4-10-II, and C3-9-R. The DFT calculations demonstrate that these four structures of OPG-Z, C4-10-I, C4-10-II, and C3-9-R are Dirac semimetals. As illustrated in Fig. 6, their linear band dispersions crossing the Fermi level at which a zero electronic density of states appears, indicating that they are the Dirac semimetals with the occurrence of massless Dirac cones. We will not discuss the OPG-Z structure because our DFT electronic band structure is in good agreement with the previous reported result.[44] For these three new Dirac semimetals, we further plot the band structure near the Dirac points in Fig. 7 to show the low-energy spectra clearer. At the Fermi level, C4-10-I, C4-10-II, and C3-9-R exhibit anisotropic Dirac cones at positions (0.021, 0.021), (0.0, 0.023), and (0.0, 0.273), respectively (the positions are in units of reciprocal vectors). In particular, for C4-10-I, the slopes of the Dirac cone are −9 eVÅ and +29 eVÅ when approaching the Dirac point from Γ to Q(0.352,0.352). For C4-10-II, the slopes are −17 eVÅ and +39 eVÅ when approaching from Γ to Y(0.0,0.5). While for C3-9-R, the slopes are 15 eVÅ and −44 eVÅ when approaching from Γ to Y(0.0,0.5). Of course, the distortion of the anisotropic Dirac cones with different slopes indicates that the electronic properties, notably conductivities, depend on directions. DFT + SOC calculations show that the SOC has no visible effects on the electronic band structures.

Fig. 7. Band structures near the Dirac points for (a) C4-10-I, (b) C4-10-II, (c) C3-9-R, and (d) C3-11. The positions of Λ are (a) (0.021, 0.021), (b) (0.0, 0.023), and (c) (0.0, 0.273). (e) and (f) The 3D visualized electronic band structures of C4-10-I and C3-9-R to highlight the shapes of the Dirac cones around the Fermi level.

3) Isotropic Dirac cones. The C3-11 structure exhibits two isotropic Dirac cones at K and K′in the hexagonal BZ (see Fig. 7(d), which is very similar to those of graphene. Their slopes are ±12 eVÅ, which is about 1/3 of that of graphene (±34 eVÅ). This means that the Dirac fermions in C3-11 have smaller Fermi velocities than those in graphene.

5. Summary

We propose two general constructing schemes of altering carbon structure on graphene. We demonstrate their usefulness by constructing numerous 2D carbon structures. We then study the structural, thermodynamic, dynamic, and electronic properties of these structures. They are all metastable compared to graphene owing to higher formation energies. Two of them ((h) OPG-A and (r) C3-10-H) are dynamically unstable, all other structures are dynamically stable. We have found one semiconductor (C3-8) and four new semimetals (C4-10-I, C4-10-II, C3-9-R, C3-11). C3-11 has isotropic Dirac cones. The other three have anisotropic Dirac cones.

In 2012, Malko et al.[27] found Dirac cones in a 2D carbon system with a rectangular unit cell. Our findings further dispute the suggestion that Dirac cones in graphene are related to its hexagonal symmetry. Indeed, out of 20 simple structures we studied, 4 new and stable structures show band crossing at the Fermi level accompanied with a vanishing electronic density. Three of them do not have a hexagonal unit cell. It shows that Dirac cones in 2D carbon systems are more common than what we previously thought.[35,48,79] Our results show that 2D carbon system is a fertile ground for searching for 2D Dirac materials.

It would be of great interest to experimentally realize these artificial two-dimensional carbon systems. The possible experimental methods include nano-scale defect-engineering on pristine graphene[36,38] and epitaxial or chemical vapor deposition on a substrate.[89,90]

Reference
[1] Novoselov K S Geim A K Morozov S V Jiang D Zhang Y Dubonos S V Grigorieva I V Firsov A A 2004 Science 306 666
[2] Novoselov K S Geim A K Morozov S V Jiang D Katsnelson M I Grigorieva I V Dubonos S V Firsov A A 2005 Nature 438 197
[3] Zhang Y Tan Y W Stormer H L Kim P 2005 Nature 438 201
[4] Neto A H C Guinea F Peres N M R Novoselov K S Geim A K 2005 Rev. Mod. Phys. 81 109
[5] Grüneis A Attaccalite C Rubio A Vyalikh D V Molodtsov S L Fink J Follath R Eberhardt W Büchner B Pichler T 2009 Phys. Rev. 80 075431
[6] Yu H L Zhai Z Y 2020 Chin. Phys. Lett. 33 117305
[7] Yin Y H Niu Y X 2020 Chin. Phys. Lett. 33 57202
[8] An Y F Dai Z H 2020 Chin. Phys. Lett. 34 17302
[9] Song L L Zhang L Z Guan Y R Lu J C Yan C X Cai J M 2019 Chin. Phys. 28 37101
[10] Zhang X F Liu Z H Liu W L Lu X L Li Z J Yu Q K Shen D W Xie X M 2019 Chin. Phys. 28 86103
[11] Zhang W Chen K B Chen Z D 2018 Acta Phys. Sin. 67 237301 in Chinese
[12] Yang Y Chen S Li X B 2018 Acta Phys. Sin. 67 237101 in Chinese
[13] Chen R K Yang C Jia Y P Gup L W Chen J N 2019 Chin. Phys. 28 117302
[14] Bolotin K I Ghahari F Shulman M D Stormer H L Kim P 2009 Nature 462 196
[15] Bonaccorso F Colombo L Yu G Stoller M Tozzini V Ferrari A C Ruoff R S Pellegrini V 2015 Science 347 1246501
[16] Zhu J Yang D Yin Z Yan Q Zhang H 2015 Small 10 3480
[17] Raccichini R Varzi A Passerini S Scrosati B 2015 Nat. Mater. 14 271
[18] Zhong S Y Shi J Luo W W Lei X L 2019 Chin. Phys. 28 78201
[19] Lin X F Zhang Z Y Yuan Z K Li J Xiao X F Hong W Chena X D Yua D S 2020 Chin. Chem. Lett. 27 1259
[20] Yoon J Sung H Lee G Cho W Ahn N Jung H S Choi M 2017 Energy Environ. Sci. 10 337
[21] Akinwande D Petrone N Hone J 2014 Nat. Commum. 5 5678
[22] Sato S 2015 Jpn. J. Appl. Phys. 54 040102
[23] Cao Y Fatemi V Demir A Fang S Tomarken S L Luo J Y Sanchez-Yamagishi J D Watanabe K Taniguchi T Efthimios K Ashoori R C Jarillo-Herrero P 2018 Nature 556 80
[24] Cao Y Fatemi V Fang S Watanabe K Taniguchi T Efthimios K Jarillo-Herrero P 2018 Nature 556 43
[25] Xu F Zhang L 2019 Chin. Phys. 28 117403
[26] Baughman R H Eckhardt H 1987 J. Chem. Phys. 87 6687
[27] Malko D Neiss C Viñes F Görling A 2012 Phys. Rev. Lett. 108 086804
[28] Kim B G Choi H J 2012 Phys. Rev. 86 115435
[29] Ivanovskii A L 2013 Prog. Solid State Chem. 4 1
[30] Narita N Nagai S Suzuki S Nakao K 1998 Phys. Rev. 58 11009
[31] Zhu H Y Balaban A T Klein D J Zivkovic T P 1994 J. Chem. Phys. 101 5281
[32] Crespi V H Benedict L X Cohen M L Louie S G 1996 Phys. Rev. 53 R13303
[33] Tyutyulkov N Dietz F Müllen K Baumgarten M 1997 Chem. Phys. Lett. 272 111
[34] Terrones H Terrones M Hernandez E Grobert N Charlier J C Ajayan P M 2000 Phys. Rev. Lett. 84 1716
[35] Wang Z Zhou X F Zhang X Zhu Q Dong H Zhao M Oganov A R 2015 Nano Lett. 15 6182
[36] Lusk M T Carr L D 2008 Phys. Rev. Lett. 100 175503
[37] Stone A J Wales D J 1986 Chem. Phys. Lett. 128 501
[38] Lusk M T Carr L D 2009 Carbon 47 2226
[39] Enyashin A N Ivanovskii A L 2011 Phys. Status Solodi 248 1879
[40] Appelhans D J Lin Z Lusk M T 2010 Phys. Rev. 82 073410
[41] Hu M Tian F Zhao Z Huang Q Xu B Wang L M Wang H T Tian Y He J 2012 Phys. Rev. Lett. 108 225505
[42] Sharma B R Manjanath A Singh A K 2014 Sci. Rep. 4 7164
[43] Lu H Li S D 2013 J. Mater. Chem. 1 3677
[44] Su C Jiang H Feng J 2013 Phys. Rev. 87 075453
[45] Wang X Q Li H D Wang J T 2013 Phys. Chem. Chem. Phys. 15 2024
[46] Wang Y Lv J Zhu L Ma Y 2010 Phys. Rev. 82 094116
[47] Glass C W Oganov A R Hansen N 2006 Comp. Phys. Comm. 175 713
[48] Xu L C Wang R Z Miao M S Wei X L Chen Y P Yan H Lau W M Liu L M Ma Y M 2014 Nanoscale 6 1113
[49] Zhang S Zhou J Wang Q Chen X Kawazoe Y Jena P 2015 Proc. Natl. Acad. Sci. USA 112 2372
[50] Ewels C P Rocquefelte X Kroto H W Rayson M J Briddon P R Heggie M I 2015 Proc. Natl. Acad. Sci. USA 112 15609
[51] Bolotin K Sikes K Jiang Z Klima M Fudenberg G Hone J Kim P Stormer H 2008 Solid State Commun 146 351
[52] Wang Z J Sun Y Chen X Q Franchini C Xu G Weng H M Dai X Fang Z 2012 Phys. Rev. 85 195320
[53] Cheng X Li R Sun Y Chen X Q Li D Li Y 2014 Phys. Rev. 89 245201
[54] Wang Z J Weng H M Wu Q S Dai X Fang Z 2013 Phys. Rev. 88 125427
[55] Li R H Ma H Cheng X Y Wang S L Li D Z Zhang Z Y Li Y Y Chen X Q 2020 Phys. Rev. Lett. 117 096401
[56] Takane D Souma S Nakayama K Nakamura T Oinuma H Hori K Horiba K Kumigashira H Kimura N Takahashi T Sato T 2018 Phys. Rev. 98 041105
[57] Li S Guo Z Fu D Pan X C Wang J Ran K Bao S Ma Z Cai Z Wang R Yu R Sun J Song F Wen J 2018 Science Bulletin 63 535
[58] Li J Ullah S Li R Liu M Cao H Li D Li Y Chen X Q 2019 Phys. Rev. 99 165110
[59] Li J Xie Q Ullah S Li R Ma H Li D Li Y Chen X Q 2018 Phys. Rev. 97 054305
[60] Li J Ma H Xie Q Feng S Ullah S Li R Dong J Li D Li Y Chen X Q 2018 Sci. China Mater. 61 23
[61] Xie Q Li J Ullah S Li R Wang L Li D Li Y Yunoki S Chen X Q 2019 Phys. Rev. 99 174306
[62] Feng B Fu B Kasamatsu S Ito S Cheng P Liu C C Feng Y Wu S Mahatha S K Sheverdyaeva P Moras P Arita M Sugino O Chiang T C Shimada K Miyamoto K Okuda T Wu K Chen L Yao Y Matsuda I 2017 Nat. Commun. 8 1007
[63] Gao L Sun J T Lu J C Li H Qian K Zhang S Zhang Y Y Qian T Ding H Lin X Du S Gao H J 2018 Adv. Mater. 30 1707055
[64] Li J Xie Q Liu J Li R Liu M Wang L Li D Li Y Chen X Q 2020 Phys. Rev. 101 024301
[65] Li J Wang L Liu J Li R Zhang Z Chen X Q 2020 Phys. Rev. 101 081403 R
[66] Ullah S Wang L Li J Li R Chen X Q 2019 Chin. Phys. 28 077105
[67] Liu Z Lou R Guo P Wang Q Sun S Li C Thirupathaiah S Fedorov A Shen D Liu K Lei H Wang S 2018 Phys. Rev. X 8 031044
[68] Chan Y H Chiu C K Chou M Y Schnyder A P 2020 Phys. Rev. 93 205132
[69] Xu Q Yu R Fang Z Dai X Weng H 2017 Phys. Rev. 95 045136
[70] Schoop L M Ali M N Straer C Topp A Varykhalov A Marchenko D Duppel V Parkin S S P Lotsch B V Ast C R 2020 Nat. Commun. 7 11696
[71] Hu J Tang Z Liu J Zhu Y Wei J Mao Z 2017 Phys. Rev. 96 045127
[72] Schilling M B Schoop L M Lotsch B V Dressel M Pronin A V 2017 Phys. Rev. Lett. 119 187401
[73] Pezzini S van Delft M R Schoop L M Lotsch B V Carrington A Katsnelson M I Hussey N E Wiedmann S 2018 Nat. Phys. 14 178
[74] Li R Li J Wang L Liu J Ma H Song H F Li D Li Y Chen X Q 2019 Phys. Rev. Lett. 123 136802
[75] Yan Z Huang P W Wang Z 2020 Phys. Rev. 93 085138
[76] Rhim J W Kim Y B 2015 Phys. Rev. 92 045126
[77] Huh Y Moon E G Kim Y B 2020 Phys. Rev. 93 035138
[78] Sankar R Peramaiyan G Muthuselvam I P Butler C J Dimitri K Neupane M Rao G N Lin M T Chou F C 2017 Sci. Rep. 7 40603
[79] Wang J Deng S Liu Z Liu Z 2015 Natl. Sci. Rev. 2 22
[80] Meille S Garboczi E J 2001 Modelling Simul. Mater. Sci. Eng. 9 371
[81] Perdew J P Burke K Ernzerhof M 1996 Phys. Rev. Lett. 77 3865
[82] Kresse G Furthmüller J 1996 Phys. Rev. 54 11169
[83] Kresse G Joubert D 1999 Phys. Rev. B, 59 1758
[84] Monkhorst H J Pack J D 1976 Phys. Rev. 13 5188
[85] Togo A Tanaka I 2015 Scr. Mater. 108 1
[86] Sobhit S Irais V J Olivia P Aldo H R 2018 Phys. Rev. 97 054108
[87] Kumar S Parks D M 2020 Journal of the Mechanics and Physics of Solids 86 19
[88] Marchenko D Varykhalov A Scholz M Bihlmayer G Rashba E Rybkin A Shikin A Rader O 2012 Nat. Commum. 3 1232
[89] Li X Cai W An J Kim S Nah J Yang D Piner R Velamakanni A Jung I Tutuc E Banerjee S K Colombo L Ruoff R S 2009 Science 324 1312
[90] Li G X Li Y L Liu H B Guo Y B Li Y J Zhu D B 2010 Chem. Commun. 46 3256
[91] Maździarz M Mrozek A Kuś W Burczyński T 2018 Materials 11 432
[92] Babu R Hiroshi M 2018 Carbon 137 266
[93] Gaikwad P V ACS Omega 4 5002
[94] Dai Z Jin W Grady M Sadowski J T Dadap J I Osgood R M Pohl K 2017 Surface Science 660 16
[95] Pop E Varshney V Roy A K 2012 MRS Bull. 37 1273